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1  Introduction 


In  this  paper  we  are  interested  in  the  Euler  equations,  the  one  dimensional  version  for  the 
perfect  gas  being  given  by 


Wt  +  f(w)x 

=  0, 

t  >  0,  x  G  R, 

(1.1) 

/  p  \ 

(  m  \ 

w  =  m  , 

f(w) 

=  pu 2  +  p 

(1.2) 

\E 

\(E  +  p)u  1 

with 

X 

m  =  pu,  E  =  — pu 2  +  pe,  p  =  (7  —  1  )pe, 

where  p  is  the  density,  u  is  the  velocity,  m  is  the  momentum,  E  is  the  total  energy,  p  is  the 
pressure,  e  is  the  internal  energy,  and  7  >  1  is  a  constant  (7  =  1.4  for  the  air).  The  speed 
of  sound  is  given  by  c  =  \fypj~p  and  the  three  eigenvalues  of  the  Jacobian  f'(w)  are  u  —  c, 
u  and  u  +  c. 

In  a  conservative  numerical  scheme,  the  internal  energy  is  obtained  by  subtracting  the 
kinetic  energy  from  the  total  energy,  thus  the  resulting  pressure  may  be  negative,  for  example, 
for  problems  in  which  the  dominant  energy  is  kinetic.  Negative  density  may  often  emerge  in 
computing  blast  waves.  Physically,  the  density  p  and  the  pressure  p  should  both  be  positive. 
The  eigenvalues  of  the  Jacobian  will  become  imaginary  if  density  or  pressure  is  negative 
so  the  initial  value  problem  for  the  linearized  system  will  be  ill-posed.  This  explains  why 
failure  of  preserving  positivity  of  density  or  pressure  may  cause  blow-ups  of  the  numerical 
algorithm. 

Replacing  the  negative  density  or  negative  pressure  by  positive  ones  is  neither  a  conser¬ 
vative  cure  nor  a  stable  solution.  Therefore,  it  is  highly  important  to  design  a  conservative 
positivity-preserving  scheme.  First  order  and  second  order  positivity-preserving  schemes 
were  well  studied,  e.g.,  [4,  11],  A  general  framework  for  constructing  arbitrarily  high  order 
positivity-preserving  discontinuous  Galerkin  (DG)  and  finite  volume  schemes  was  proposed 
recently  in  [20].  This  framework  can  be  easily  generalized,  for  instance,  to  unstructured 
meshes  [23],  and  to  general  equations  of  state  and  Euler  system  with  source  terms  [21]. 
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Generalization  of  the  positivity-preserving  method  in  [20]  to  high  order  finite  difference 
schemes  is  not  straightforward.  On  the  other  hand,  in  some  applications  where  high  or¬ 
der  schemes  are  preferred,  for  example,  cosmological  simulation  [5],  finite  difference  WENO 
scheme  [10]  is  more  favored  than  DG  schemes  [2,  3]  and  the  finite  volume  WENO  scheme 
[12,  15]  due  to  its  smaller  memory  cost  (compared  to  DG)  and  smaller  computational  cost 
(compared  both  to  finite  volume  schemes  and  to  DG  schemes)  for  multi-dimensional  prob¬ 
lems,  see  for  example  a  comparison  in  [1]  in  the  context  of  ENO  schemes. 

In  this  paper,  we  will  follow  the  idea  in  [20]  to  construct  positivity-preserving  high  order 
finite  difference  WENO  schemes.  We  will  show  that  by  adopting  the  same  simple  limiter  as 
in  [20],  to  a  slightly  different  version  of  finite  difference  WENO  schemes  from  the  one  in  [10], 
the  final  scheme  will  keep  the  positivity  of  density  and  pressure  without  losing  conservation. 

A  conservative  positivity-preserving  scheme  is  L1-stable,  see  [22],  The  limiter  will  not  destroy 
the  high  order  accuracy  of  the  WENO  scheme  for  smooth  solutions  without  vacuum.  All  the 
results  also  hold  for  finite  difference  ENO  schemes  [16]. 

The  paper  is  organized  as  follows.  In  Section  2,  we  briefly  review  the  positivity-preserving 
finite  volume  schemes  in  [20]  and  the  finite  difference  WENO  scheme  in  [10].  Then  we 
introduce  positivity-preserving  finite  difference  schemes  in  one  space  dimension  for  the  perfect 
gas  in  Section  3.  In  Section  4,  we  discuss  a  straightforward  extension  to  multi-dimensions, 
general  equations  of  state  and  source  terms.  In  Section  5,  numerical  tests  of  the  fifth  order 
WENO  schemes  for  some  very  demanding  problems  are  shown.  Concluding  remarks  are 
given  in  Section  6. 

2  Preliminaries 

2.1  Review  of  positivity-preserving  high  order  finite  volume  WENO 
schemes 

We  first  briefly  review  the  basic  idea  in  [20,  22]  for  finite  volume  WENO  schemes.  Consider 
the  Euler  equations  (1.1)  in  more  detail.  Let  p( w)  =  (7  —  1  )(E  —  be  the  pressure 
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function.  It  can  be  easily  verified  that  p  is  a  concave  function  of  w  =  ( p,m,E)T  if  p  >  0. 
For  wi  =  (pi,rrii,  E1)T  and  w2  =  (p2,  nr2,  E2)T,  Jensen’s  inequality  implies,  for  0  <  s  <  1, 

p  (swi  +  (1  —  s)w2)  >  sp  (wi)  +  (1  —  s)p  (w2) ,  if  Pi  >  0,  p2  >  0.  (2.1) 


Define  the  set  of  admissible  states  by 


G  = 


p  >  0  and  p 


(7  -  1) 


1  m2\ 

2~) 
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then  G  is  a  convex  set.  We  want  to  construct  finite  volume  WENO  schemes  producing 
solutions  in  the  set  G. 

The  time  discretization  is  chosen  as  the  high  order  strong  stability  preserving  (SSP) 
methods  [14,  16,  8,  9]  which  are  convex  combinations  of  Euler  forward.  Thus  we  only  need 
to  discuss  the  Euler  forward  since  G  is  convex. 

A  general  high  order  finite  volume  scheme  has  the  following  form 


w-  =  w, 


f  I  W  !  ,  W 


w 


2  1  2 


(2.2) 


where  f  is  a  positivity  preserving  flux,  for  instance,  Lax-Friedrichs  flux,  w”  is  the  approxi¬ 
mation  to  the  cell  average  of  the  exact  solution  v(x,t)  in  the  cell  =  [Xj_i,o;i+i]  at  time 
level  n,  and  w“  i,  w^i  are  the  high  order  approximations  of  the  point  values  v(xi+i ,  tn) 
within  the  cells  and  It+\  respectively.  These  values  are  reconstructed  from  the  cell  av¬ 
erages  w"  by  the  WENO  reconstruction.  We  assume  that  there  is  a  polynomial  vector 
qi(x)  =  (pi(x),  rni(x),  Ei(x))T  with  degree  k  which  are  (k  +  l)-th  order  accurate  approxi¬ 
mations  to  smooth  exact  solutions  v(x,t)  on  A,  and  satisfies  that  w”  is  the  cell  average 
of  qj(a;)  on  Ii:  =  q^x^i)  and  wT,  =  qj(xi+i).  The  existence  of  such  polynomials 

can  be  established  by  interpolation  for  WENO  schemes.  For  example,  for  the  fifth  or¬ 
der  WENO  scheme,  there  is  a  unique  vector  of  polynomials  of  degree  four  q,;(x)  satisfying 
q«  0  i- 1 )  =  w+  1 ,  q*  {xi+ 1 )  =  w~+ 1  and 

f  q i(x)dx  =  w j,  for  j  —  i  —  1,  i,  i  4- 1. 

Ii 


1 

Ax 
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We  need  the  Appoint  Legendre  Gauss-Lobatto  quadrature  rule  on  the  interval  = 
[%i_ i,xi+i],  which  is  exact  for  the  integral  of  polynomials  of  degree  up  to  2N  —  3.  We 
would  need  to  choose  N  such  that  2N  —  3  >  k.  Denote  these  quadrature  points  on  Ll  as 


Si  =  {aq  i  =  x],x2,---  ,x\ 


N~l  x? 


(2.3) 


In  [20],  we  have  shown  that  q i(xf)  G  G  for  all  j  and  a  is  a  sufficient  condition  for 
w"+1  G  G  under  suitable  CFL  conditions.  The  limiter  in  [20]  can  enforce  this  sufficient 
condition  without  destroying  conservation.  Moreover,  the  limiter  is  an  accurate  modification 
for  smooth  solutions  if  there  is  no  vacuum  region  in  the  exact  solution.  See  [22,  17]  for  simpler 
implementations  of  the  limiter. 


2.2  Review  of  the  finite  difference  WENO  scheme  for  scalar  con¬ 
servation  laws 


Before  describing  the  finite  difference  WENO  scheme  for  the  Euler  system,  we  first  review 
the  relations  between  finite  difference  and  finite  volume  WENO  schemes  in  [15]  for  the  scalar 
linear  equation 


ut  +  ux  =  0,  u(x,  0)  =  u0(x).  (2.4) 

Consider  a  uniform  mesh  with  nodes  ay.  Define  xi+i  =  |(ay+i  +  ay),  A  =  [xj_i,xi+i]  and 
Ax  =  Xi+ 1  —  x%.  A  conservative  finite  difference  scheme  with  high  order  spatial  discretization 
and  Euler  forward  time  discretization  solving  (2.4)  has  the  form 


u?+1  =  u? 


Arc  ^+5 


where  ^[/i+i  —  i]  should  be  a  high  order  approximation  to  ux  at  x  =  Xi. 
If  there  exists  a  function  h(x)  depending  on  the  mesh  size  Ax  such  that 


1 


i  Ax 

rx+^r 


u (x)  =  -t—  /  h(£)d£, 

Ax  ,/™_A x 


(2.5) 


(2.6) 


then  we  call  u  and  h  a  reconstruction  pair.  See  [6]  for  a  detailed  discussion  of  the  recon¬ 
struction  pair.  Denote  them  by 


h  =  Rax(u),  u  = 
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Then  ux  =  [h  (x  +  —  h  [x  —  4r)]  .  Thus  if  fi+i  is  a  high  order  accurate  approxima¬ 

tion  of  h(xi+ 1)  then  -^[fi+ 1  —  /)_i]  will  be  a  high  order  approximation  to  ter  at  x  =  Xi. 

Notice  that  the  cell  average  of  h  over  I,  is  simply  the  point  value  u(xi)  by  (2.6).  So  the 
high  order  finite  difference  WENO  scheme  for  (2.4)  can  be  formulated  as 

•  At  time  level  n,  obtain  the  cell  averages  of  h  on  /*  by  1\  =  u”  . 


•  Use  the  fifth  order  WENO  reconstruction  based  on  h"  (j  in  a  neighborhood  of  i)  to 
construct  nodal  values  of  h  at  x; ,  1 ,  and  denote  it  by  h~  1 . 

*+2  »+5 

•  Set  the  flux  as  fi+ 1  =  /iT^i  in  (2.5). 


It  is  clear  that  (2.5)  is  equivalent  to 


-r-n+1  — n 

ru  =  h. 


Aa: 


h: 


.1  j> 


(2.7) 


which  is  a  finite  volume  scheme  for  the  function  h.  In  other  words,  for  (2.4),  a  finite  difference 
WENO  scheme  of  u  is  exactly  a  finite  volume  WENO  scheme  of  h. 

Let  qi{x)  denote  a  polynomial  which  is  a  high  order  accurate  approximation  to  h  on  J* 
and  satisfies  qi(xi+ 1)  =  h~+l.  By  the  result  in  [19],  if  qi(x f)  G  [m,  M]  where  m  and  M  are 
the  minimum  and  maximum  of  the  initial  value  Uq(x),  then  hi  G  [m,  M]  under  certain 
CFL  conditions.  It  appears  that  we  can  easily  construct  maximum-principle-satisfying  finite 
difference  WENO  schemes  for  scalar  conservation  laws  following  [19].  Unfortunately,  this 
method  will  destroy  accuracy.  Even  though  h  will  converge  to  u  when  Ax  goes  to  zero,  for 
a  fixed  Ax,  h  has  larger  maximum  and  smaller  minimum  than  u  does  by  the  analysis  in  [6] . 
Since  qt{x)  approximates  h  rather  than  u,  enforcing  q.i(x f)  G  [m,  Ad]  will  destroy  the  high 
order  accuracy. 

However,  if  the  minimum  of  u(x,t)  is  strictly  positive,  then  h(x,t)  >  0  for  small  enough 
Ax,  thus  enforcing  qt(xf)  >  0  will  not  destroy  high  order  accuracy.  Therefore,  we  can  still 
take  advantage  of  this  relation  to  construct  high  order  positivity-preserving  finite  difference 
WENO  scheme  in  the  next  section. 
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3  Positivity-preserving  high  order  finite  difference  WENO 
schemes  in  one  dimension 

3.1  A  finite  difference  WENO  scheme  for  Euler  equations 


We  now  discuss  a  (k  +  l)-th  order  finite  difference  WENO  scheme  for  (1.1).  Given  the 
point  values  w™  at  time  level  n.  let  Ai+i  denote  the  Roe  matrix  [13]  of  the  two  states  w”+1 
and  w”.  Let  Li+ 1  and  Ri+i  be  the  left  and  right  eigenvector  matrices  of  Ai+i  respectively, 

i.e.,  A  =  RAL  where  A  is  the  diagonal  matrix  with  eigenvalues  of  A  on  the  diagonal.  Let 
a  =  max  ||  (|tt|  +  c)  ||  where  u  and  c  are  the  velocity  and  speed  of  sound  of  the  state  w"  and 
the  maximum  is  taken  either  globally  over  all  the  w"  or  locally  over  the  w"  in  the  WENO 
reconstruction  stencil.  We  will  use  the  Lax- Friedrichs  splitting, 


(3.1) 


At  each  fixed  xi+i, 

1.  Let  h±  =  -Ra^P)-  Then  we  have  the  cell  averages  h±”'  =  f±(w)1). 

2.  Transform  all  the  cell  averages  h±”  (J  in  a  neighborhood  of  i)  to  the  local  characteristic 
fields  by  setting 


_  Yi  T  i  71 

U±j  =  • 

3.  Perform  the  WENO  reconstruction  for  each  component  of  u+”  to  obtain  approxima¬ 
tions  of  the  point  value  of  the  function  Li+  ih+  at  the  point  xi+i  and  denote  them 
as  (u+)^j_.  Perform  the  WENO  reconstruction  for  each  component  of  u_"  to  obtain 
approximations  of  the  point  value  of  the  function  Li+ 1  h_  at  the  point  xi+i  and  denote 
them  as  (u_)^_!.  Here  the  superscripts  —  and  +  denote  approximations  within  the 
cells  Ii  and  Il+\  respectively. 

4.  Transform  back  into  physical  space  by 
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(3.2) 


(h+) 


Form  the  flux  by 


OL 


(hn 


>i+l 


Let  A  =  At/ Ax,  then  we  get  the  conservative  scheme 


w”+‘ =  W?  -  A  (t+.  -?_>).  (3.3) 

3.2  A  sufficient  condition  to  keep  positivity 

Next,  we  will  derive  a  sufficient  condition  for  the  scheme  (3.3)  to  keep  w”+1  e  G  provided 
w”  G  G.  (3.1)  and  (3.2)  imply  that  (3.3)  can  be  written  as 


w 


,71+1 


w»-A(f)+.-U.) 

h+1  +  -  Aa  ((h+)“+,  -  (h_)+  ,  -  (h+)‘_ ,  +  (h_)+  .) 

H+  +  H-, 


where 


H+ 


aX 
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H 


h,  +  aX 


(h-) 


(3.4) 

(3.5) 


Notice  that  (3.4)  and  (3.5)  are  two  finite  volume  WENO  schemes  for  h+  and  h_.  Obvi¬ 
ously  it  suffices  to  discuss  conditions  to  keep  H±  e  G.  We  first  discuss  (3.4).  By  interpola¬ 
tion,  there  exists  a  vector  of  polynomials  of  degree  k  qf(x)  such  that  q+(xi+i)  =  (h+)7^i, 
the  cell  average  of  q+(xi+i)  on  is  h+”  and  q+(x)  is  a  (k  +  l)-th  order  accurate  approx¬ 
imation  to  the  function  h+  on  the  cell  if  w  is  smooth.  The  exactness  of  the  quadrature 
rule  for  polynomials  of  degree  k  implies 


N 

q  t(x)dx  =  ^2waqt(xf) 

a=l 


(3.6) 


N- 1 

=  $<*<£&)  +  Wiv(h+)-+i  =  (1  -  wN)qp*  +  ^(h+)“+i, 

a=l 

N-l 

where  we  use  q^'+  =  1_~L  Y)  waq f(xf)  to  denote  a  convex  combination  of  point  values  of 

a=l 

Plugging  (3.6)  into  (3.4),  we  have 


H+  =  (1  -  wN)qp*  +  ( wN  -  «A)(h+).+  i  +  aA(h+)._i. 

Therefore,  if  q]4'*,  (h+)7  ±,  (h+)~_i  G  G,  then  H+  G  G  under  the  CFL  condition  aX  <  wn- 

l+2  1  2 

Similarly,  for  (3.5),  there  exists  a  vector  of  polynomials  of  degree  k  q“(a;)  such  that 

qi(x^i)  =  (h_)^i,  the  cell  average  of  q^(Xj_i)  on  I,  is  and  q~(a;)  is  a  (k  +  1)- 
2*2  2 

th  order  accurate  approximation  to  the  function  h_  on  the  cell  /,  if  w  is  smooth.  Let 

N 

qr*  =  tzW  E  Wa<h(xf),  then  (3.5)  becomes 

a=2 

H"  =  (1  -  Wi)q~’*  +  (w1  -  aA)(h_)+ i  +  aA(h_)+  . 

1  2  Z_h2 

Therefore,  if  q,*,  (h_)+  i,  (h_)A_i  G  G,  then  H  G  G  under  the  CFL  condition  aX  <  vd\. 

*+2  *  2 

Notice  that  w i  =  wn,  we  have 


Theorem  3.1  Under  the  CFL  condition  aX  <  W\,  if  q4-’*,  (h+) .  i;  (h+)  i,qi  (h.)4,  k, 

1  '  2  1  2  1  '  2 

(h_)j^i  G  G,  the  finite  difference  WENO  scheme  (3.2)  and  (3.3)  will  be  positivity-preserving, 
1  2 

i.e.,  w”+1  G  G,  where  by  (3.6)  we  have 


q 


+,* 


l 

1  -  wn 


«hv(h+)i+ij  , 


1 

1  —  W\ 


wi(h_)+i 
1  2 


(3.7) 


3.3  A  simple  limiter 


At  time  level  n,  given  w"  G  G  for  all  i,  then  h±™  =  f±(w”)  G  G,  which  was  proved  in 
Remark  2.4  of  [20].  Thus  we  can  enforce  the  positivity  of  the  two  finite  volume  schemes 
(3.4)  and  (3.5)  by  the  techniques  for  finite  volume  WENO  schemes  in  [20,  22], 

For  convenience,  we  use  p{ w)  and  p( w)  to  denote  the  density  and  pressure  of  a  state 


w.  Let  h+-  = 


=  (pi,m;,£A)T  (h+).+  i  =  (pi+i,m.+  i,E.+^  and  q+’*  =  ( p*,m*,E*)r 
explicit  formula  of  q f(x)  will  not  be  needed.  Consider  the  following  limiter  for  (3.4), 


The 
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Setup  a  small  number  £  =  min  jlO  13,  p  }• 

For  each  cell  A,  modify  the  density  first: 

i  =  0i  fpZ-i  “  Pi)  +  Pn  9i=  min 


Pi  ~  g 
Pi  Pmir 


■a  , 


(3.8) 


where  pmin  =  min  | Pi+i  ,  P* )  •  Notice  that  both  the  numerator  and  the  denominator 
of  are  nonnegative.  Then  denote  (h+)T,  \  —  \p~,  i,m7,  i ,E7,  i  )  and  q.+’*  = 

Pi  Pmin  v  7  7--* —  '  '>■-* —  ‘  ,  —  •  L 


1 


l-wN 


h+r  -  Wiv(h+)  i 


Then  modify  the  pressure  for  each  cell  For  convenience,  let  q1  =  (h+).+i  and 
q2  =  q^1"’*.  For  m  =  1,  2,  if  p( qm)  <  0,  then  solve  the  following  equation  for  tm  G  [0, 1), 


P 


(1  -  tm) h+”  +  tmqn 


=  0. 


(3.9) 


Notice  that  the  convexity  of  G  ensures  the  uniqueness  of  tm  G  [0,1).  Although  (3.9) 
is  only  a  quadratic  equation  for  the  ideal  gas,  it  is  more  robust  to  avoid  solving  it  due 
to  round  off  errors.  By  the  discussion  in  Section  3  of  [17],  instead  we  can  solve  the 
following  linear  equation  for  tm, 


(1  -  tm)p  (h+” )  +  tmp  (qm)  =  0. 


(3.10) 


Then  Jensen’s  inequality  (2.1)  implies  that  the  solution  of  (3.10)  satisfies 


P 


(1  -  tm) h+”  +  tm q” 


>  0. 


If  p( qm)  >  0,  set  tm  =  1.  Get  the  modified  point  value  by 


(h+)i+i  =  ( (h+)  i  -  h+”)  +  h+”,  d2  =  min  {t\  t2}  . 


(3.11) 


The  limiter  for  (3.5)  can  be  defined  in  the  same  way.  Similarly,  we  get  the  revised  point 

value  (h_)j^i.  Then  we  have  the  modified  WENO  scheme  (3.3)  with 

l~2 


(3.12) 
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It  is  straightforward  to  check  that  the  new  scheme  (3.3)  and  (3.12)  satisfies  the  sufficient 
condition  in  Theorem  3.1. 

Suppose  v(x,t)  is  a  smooth  solution  of  (1.1)  and  it  satisfies  minp(v(a:,  t))  >  0  and 

X,t 

minp(v(:r, t))  >  0,  then  (v(z,t))  =  JRA.r(f±(v(  x,t)))  G  G  for  any  (x,t)  if  Ax  is  small 

X,t 

enough.  Since  the  limiter  (3.8)  and  (3.11)  is  the  exactly  the  same  limiter  for  finite  volume 
scheme  (3.4)  as  in  [20,  22],  following  the  same  arguments  in  [20],  it  is  straightforward  to 
show  that  the  accuracy  will  not  be  destroyed  by  the  limiter  for  smooth  solutions  without 
vacuum  regions  when  Ax  is  small. 

Remark  3.2  The  limiter  should  be  used  for  each  stage  in  a  SSP  Runge-Kutta  method  or  each 
step  in  a  SSP  multi-step  method.  To  save  computational  cost,  the  stringent  CFL  condition 
aX  <  W\  can  be  strictly  enforced  only  when  a  precalculation  to  the  next  time  stage  or  time 
step  produces  negative  density  or  negative  pressure. 


4  Generalizations 

4.1  Source  terms 


The  one  dimensional  version  of  the  Euler  equations  with  source  terms  is  given  by 


wf  +  f(w)x  =  s(w,  x),  t  >  0,  x  G  It, 
The  finite  difference  WENO  scheme  for  (4.1)  can  be  written  as 


w”+1  =  wj 


a  I Gi 


fi-i  )  +  Afs(w?n,  Xi)  =  +  is, 


with 


I'  _  w"  -  2A  f- ,  i  -  £ 


(4.1) 


(4.2) 

S  =  wf  +  2Ats(wf,a;i).  (4.3) 

For  positivity  of  w”+1,  it  suffices  to  have  F.SgG.  (4.2)  is  in  the  same  form  as  (3.3),  thus 
the  positivity  of  F  is  straightforward.  It  was  discussed  in  [21]  that,  for  four  kinds  of  source 
terms  (geometric,  gravity,  chemical  reaction  and  radiative  cooling),  S  G  G  under  a  suitable 
CFL  condition  if  w”  G  G. 
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4.2  General  equations  of  state 

Consider  a  general  equation  of  state  satisfying  the  following  assumption: 


if  p  >  0,  then  e  >  0  p  >  0.  (4.4) 

_ 1  rr? 

The  pressure  is  not  necessarily  a  concave  function.  However,  the  internal  energy  e  =  — ^-e- 
is  always  a  concave  function  of  w,  thus  G  =  {w  :  p  >  0,p  >  0}  is  still  convex. 

Therefore,  it  is  straightforward  to  check  that  Theorem  3.1  still  holds.  By  setting  a  > 
max  ||  |u|  +  ||  ,  it  is  still  true  that  f^w)  G  G  if  w  G  G,  see  [21].  Replace  (3.9)  by 


e 


(1  -  tm) h+;  +  tmq 


0, 


or  replace  (3.10)  by 

(1  -  tm)e  (h+")  +  tme  (qm)  =  0, 

then  the  same  limiter  (3.8)  and  (3.11)  can  be  used  to  enforce  the  sufficient  condition. 


4.3  Extensions  to  multi-space  dimensions 


Two  dimensional  Euler  equations  are  given  by 


with 


W,  +  f(w),  +  g(w),  =  0,  t>  0,(i,i)eB’ 


( p  \ 

/  m  \ 

/  n  ^ 

m 

,  f(w)  = 

pu 2  +  p 

,  g(w)  = 

puv 

9  1 

n 

puv 

pv  +  p 

E  j 

\  ( E  +  p)u  ) 

\(E+  p)v  / 

m  —  pu ,  n  —  pv,  E  —  +  Pe>  P  —  (7  —  l)Pe> 


(4.5) 


(4.6) 


where  p  is  the  density,  u  is  the  velocity  in  x  direction,  v  is  the  velocity  in  y  direction,  m  and 
n  are  the  momenta,  E  is  the  total  energy,  p  is  the  pressure,  e  is  the  internal  energy.  The 
speed  of  sound  is  given  by  c  =  \fypjp-  The  eigenvalues  of  the  Jacobian  f;(w)  are  u  —  c,  u,  u 
and  u  +  c  and  the  eigenvalues  of  the  Jacobian  g'(w)  are  v  —  c,  v,  v  and  v  +  c.  The  pressure 
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G  =  <  w  = 


P>  0  and  p=  (7^  1)  (e-  >  0 


function  p  is  still  concave  with  respect  to  w  if  p  >  0  and  the  set  of  admissible  states 

(  P  \ 

m 
n 

\  E  J 

is  still  convex. 

We  will  use  the  same  Lax-Friedrichs  splitting 

,  ay  =  max  ||  (|u|  +  c)  ||, 

,  a2  =  max  ||  (|u|  +  c)  ||, 

where  the  maximum  is  taken  either  globally  or  locally.  Following  Remark  2.4  in  [20],  it  is 
easy  to  check  that  f^w),  g±(w)  e  G  if  w  e  G. 

Assume  a  uniform  mesh  with  nodes  (xt,  ijj),  the  finite  difference  WENO  scheme  is  given 
by 

=  (t+U  -<W)  -  |£  (g.J+i  -iij-i) . 

where  the  fluxes  f \+i  j  and  g,  ]+±  are  obtained  by  the  same  one- dimensional  WENO  approx¬ 
imation  as  in  the  previous  section. 

For  the  discussion  of  positivity,  rewrite  the  scheme  as  1  =  |F  +  |G  with 


F^w) 

1 

fw±f(w)' 

_  2  1 

V  “i 

g±(w) 

1 

(w±S(w) 

_  2 

V  «2 

'  w‘  j  ' 

At 


(4.7) 

G  =  Wu  -  -  <4’8) 

If  F,  G  G  G,  then  w”^1  G  G.  Notice  that  (4.7)  and  (4.8)  share  the  same  abstract  form  as 
(3.3),  it  is  straightforward  to  extend  the  positivity-preserving  results  for  one-dimension  to 
multi-space  dimensions. 

5  Numerical  tests 


In  this  section,  we  show  some  results  of  the  fifth  order  finite  difference  WENO  scheme  and 
the  third  order  Runge-Kutta  with  the  positivity-preserving  limiter  for  several  demanding 
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examples.  The  WENO  schemes  without  the  positivity-preserving  limiter  will  blow  up  for 
most  examples  in  this  section. 


Example  5.1  Accuracy  Test. 


Consider  the  vortex  evolution  problem  for  (4.5).  The  mean  flow  is  p  =  p  =  u  =  v  =  l.  Add 
to  the  mean  flow  an  isentropic  vortex  perturbation  centered  at  (x0,  yo)  in  (u,  v)  and  T  =  p/ p, 
no  perturbation  in  S  =  p/p^, 

(Su.Sv)  =  ±e^H-y,x),  ST  =  , 

where  (x,  y)  =  (x  —  x0,  y  —  yo),r2  =  x2  +  y2.  The  exact  solution  is  the  passive  convection  of 
the  vortex  with  the  mean  velocity. 

The  domain  is  taken  as  [0, 10]  x  [0, 10]  and  ( Xo,yo )  =  (5,5).  The  boundary  condition  is 
periodic.  We  set  7  =  1.4  and  the  vortex  strength  e  =  10.0828  such  that  the  lowest  density 
and  lowest  pressure  of  the  exact  solution  are  7.8  x  10~15  and  1.7  x  10”20.  We  test  the  accuracy 
of  the  limiter  (3.8)  and  (3.11)  on  the  fifth  order  finite  difference  WENO  scheme  with  the 
third  order  SSP  Runge-Kutta.  In  order  to  make  the  error  in  spatial  discretizations  dominant, 
we  take  At  =  Ax^ .  See  Table  5.1.  We  clearly  observe  the  fifth  order  accuracy. 

Table  5.1:  Example  5.1:  Fifth  order  finite  difference  WENO  scheme  with  the  positivity¬ 
preserving  limiter,  for  the  vortex  evolution  problem,  T  =  0.01,  and  Ax  =  Ay. 


1/  Ax 

L 1  error 

order 

L°°  error 

order 

8 

6.77e-6 

- 

5.33e-4 

- 

16 

3.26e-7 

4.37 

3.77e-5 

3.82 

32 

8.04e-9 

5.34 

1. Ole-6 

5.22 

64 

1.92e-10 

5.38 

3.62e-8 

4.80 

Example  5.2  ID  low  density  and  low  pressure  problems. 

We  consider  two  one-dimensional  low  density  and  low  pressure  problems  of  (1.1)  for  ideal 
gas.  The  first  one  is  a  one- dimensional  Riemann  problem,  for  which  the  initial  condition 
is  Pl  =  Pr  =  7,  ul  =  —1,  Ur  —  1,  pl  —  Pr  —  0.2  and  7  =  1.4.  The  exact  solution 
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contains  vacuum.  It  is  the  same  double  rarefaction  problem  in  [11],  The  results  of  positivity¬ 
preserving  fifth  order  finite  difference  WENO  scheme  are  shown  in  Figure  5.1  (left),  which 
are  comparable  to  the  results  of  positivity-preserving  DG  method  [20]. 

The  second  one  is  one- dimensional  Sedov  blast  wave.  For  the  initial  condition,  the  density 
is  1,  velocity  is  zero,  total  energy  is  10~12  everywhere  except  that  the  energy  in  the  center 
cell  is  the  constant  ^  with  E0  =  3200000  (emulating  a  5-function  at  the  center).  7  =  1.4. 
The  results  of  positivity-preserving  fifth  order  finite  difference  WENO  scheme  are  shown  in 
Figure  5.1  (right),  which  are  again  comparable  to  the  results  of  positivity-preserving  DG 
method  [20]. 

Example  5.3  2D  low  density  and  low  pressure  problems. 

We  consider  two  two-dimensional  low  density  and  low  pressure  problems  of  (4.5)  for  ideal  gas. 
The  Erst  one  is  two-dimensional  Sedov  blast  wave.  The  computational  domain  is  a  square. 
For  the  initial  condition,  the  density  is  1,  velocity  is  zero,  total  energy  is  10-12  everywhere 
except  that  the  energy  in  the  lower  left  corner  cell  is  the  constant  •  7  =  1.4.  The 

numerical  boundary  treatment  is  reflective  for  the  left  and  bottom  edges.  See  Figure  5.2. 

The  second  one  is  shock  diffraction  problem.  The  setup  is  the  following:  the  computa¬ 
tional  domain  is  the  union  of  [0, 1]  x  [6, 11]  and  [1, 13]  x  [0, 11];  the  initial  condition  is  a  pure 
right-moving  shock  of  Mach  =  5.09,  initially  located  at  x  —  0.5  and  6  <  y  <  11,  moving  into 
undisturbed  air  ahead  of  the  shock  with  a  density  of  1.4  and  pressure  of  1.  The  boundary 
conditions  are  inflow  at  x  —  0,  6  <  y  <  11,  outflow  at  x  —  13,  0  <  y  <  11,  1  <  x  <  13,  y  —  0 
and  0  <  x  <  13,  y  —  11,  and  reflective  at  the  walls  0  <  x  <  1,  y  —  6  and  x  —  1,  0  <  y  <  6. 
7  =  1.4.  The  density  and  pressure  at  t  =  2.3  are  presented  in  Figures  5.3. 

The  results  are  comparable  to  those  of  positivity-preserving  DG  method  [20]  and  finite 
volume  WENO  scheme  [22], 

Example  5.4  High  Mach  number  flows. 
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6 


(a)  density 


(b)  density 


(c)  pressure 


(d)  pressure 


Figure  5.1:  Example  5.2.  Left:  Double  rarefaction  problem.  Right:  ID  Sedov  blast.  Ax  = 
1/200.  The  red  curves  are  exact  solutions.  Green  symbols  are  numerical  solutions  of  fifth 
order  WENO  with  the  positivity- preserving  limiter. 
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(a)  Color  contour  of  density  (b)  Cut  along  y  =  0.  The  solid  line  is  the 

exact  solution.  Symbols  are  numerical  solu¬ 
tions. 


Figure  5.2:  2D  Sedov  blast.  T  =  1.  Ax  =  Ay  =  ^1.  The  fifth  order  finite  difference  WENO 
scheme  with  positivity-preserving  limiter. 


(a)  Density:  20  equally  spaced  contour  lines  (b)  Pressure:  40  equally  spaced  contour 
from  p  =  0.066227  to  p  =  7.0668.  lines  from  p  =  0.091  to  p  =  37. 

Figure  5.3:  Shock  diffraction  problem.  Ax  =  Ay  =  1/80.  The  fifth  order  finite  difference 
WENO  scheme  with  positivity-preserving  limiter. 
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We  consider  the  Mach  2000  problem  in  [20].  The  equations  are  (4.5)  with  7  =  1.2.  The  do¬ 
main  is  [0, 1]  x  [—0.25,  0.25],  initially  full  of  the  ambient  gas  with  (p,  u,  v,p)  =  (5,  0,  0,  0.4127). 
The  boundary  conditions  for  the  right,  top  and  bottom  are  outflow.  For  the  left  boundary, 
(p,  u,  v,  p)  =  (5,  800,  0,  0.4127)  if  y  e  [—0.05,  0.05]  and  ( p ,  u,  v,  p)  =  (5,  0,  0,  0.4127)  otherwise. 
The  terminal  time  is  0.001.  The  speed  of  the  jet  is  800,  which  is  around  Mach  2100  with 
respect  to  the  soundspeed  in  the  jet  gas.  See  Figure  5.4  for  the  result  on  a  800  x  400  mesh, 
which  is  comparable  to  the  result  in  [20]. 


Figure  5.4:  Example  5.4.  Simulation  of  Mach  2000  jet.  Log  plot  of  density. 


Example  5.5  The  reactive  Euler  equations. 


We  consider  the  following  equations  which  are  often  used  to  model  the  detonation  waves: 


w,  +  f(w)„  +  g(w)„  =  s(w),  l>0,(i,()eF’,  (5.1) 


p  \ 

/  m  \ 

(  n  ^ 

f 0  \ 

m 

pu 2  +  p 

puv 

0 

n 

,  f(w)  = 

puv 

,  g(w)  = 

pv 2  +  p 

,  s(w)  = 

0 

E 

(. E  +  p)u 

(E  +  p)v 

0 

pY  yJ 

\  puY  y 

\  PvY  / 

\w/ 

(5.2) 


with 


m  =  pu, 


n  =  pv,  E 


1  9  1  2 

-pu  +-pv  + 


P 

7-1 


+  pqY, 
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where  q  is  the  heat  release  of  reaction,  7  is  the  specific  heat  ratio  and  Y  denotes  the  reactant 
mass  fraction.  The  source  term  is  assumed  to  be  in  an  Arrhenius  form 

cu  =  -KpYe~f/T, 

where  T  =  -p  is  the  temperature,  T  is  the  activation  temperature  and  K  is  a  constant.  The 
eigenvalues  of  the  Jacobian  f'(w)  are  u  —  c,u,u,u,u  +  c  and  the  eigenvalues  of  the  Jacobian 
g'(w)  are  v  —  c,  v,  v,  v,  v  +  c,  where  c  = 

Consider  the  following  problem  which  is  similar  to  the  shock  diffraction  problem.  The 
computational  domain  is  the  union  of  [0, 1]  x  [2,  5]  and  [1,  5]  x  [0, 5];  The  initial  conditions  are, 
if  x  <  0.5,  (p,  u,  v,  E,  Y)  =  (11,  6.18,  0,  970, 1);  otherwise,  (p,  u,  v,  E,  Y)  =  (1,  0,  0,  55, 1).  The 
boundary  conditions  are  reflective  except  that  at  x  —  0,  (p,  u,  v,  E,  Y)  =  (11,  6.18,  0,  970, 1). 
The  terminal  time  is  t  —  0.6.  The  parameters  are  7  =  1.2,  q  =  50,  T  =  50,  K  =  2566.4. 

See  Figure  5.5,  which  is  comparable  to  the  result  of  the  third  order  positivity-preserving 
DG  method  in  [17]. 

Example  5.6  A  general  equation  of  state. 

We  consider  the  three  species  model  with  a  more  general  equation  of  state  in  [18].  The  model 
involves  three  species,  02,  O  and  N2  (pi  =  po,  P2  =  Po2  and  Pz  =  Pn2)  with  the  reaction: 

o2  +  n2^±o  +  o  +  n2. 


The  governing  equations  are 


(  pl\ 

(  Plu  \ 

/  2 Mxu  \ 

P2 

P'2U 

—M2oj 

P3 

+ 

P3U 
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0 

pu 

pu 2  +  p 

0 

E  ) 

t 

\  ( E  +  p)u  ) 

X 

V  0  / 

and 


p  =  5>, 

5=1 


p  = 


S=1 


Ps 

M, 


E 


^pses(T)  +  p1h'i  +  - pu 2 

S=  1 
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(a)  Colored  contour  of  density 


(b)  Contour  line  of  density 


(c)  Colored  contour  of  pressure 


Figure  5.5:  Example  5.5.  Detonation  diffraction  at  a  90°  corner.  Ax  =  Ay  =  1/80.  Fifth 
order  WENO  scheme  with  positivity-preserving  limiter. 
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where  the  enthalpy  h\  is  a  constant,  R  is  the  universal  gas  constant,  Ms  is  the  molar  mass  of 
species  s,  and  the  internal  energy  es(T)  =  3R/2Ads  and  5R/2MS  for  monatomic  and  diatomic 
species  respectively.  The  rate  of  the  chemical  reaction  is  given  by 


uj  = 


k^W2  -  k-(T)  (it )  1  tw-  kf  =  cr-U-'r 


s=l 


kb  =  kf/  exp  ( bi  +  l°g  z  +  b3z  +  b4z2  +  b5z3),  z  =  10000/T 

where  bi,  C  and  E  are  constants  which  can  be  found  in  [18,  7]. 

The  eigenvalues  of  the  Jacobian  f'(w)  are  {u,u,u,u  +  c,u  —  c)  where  c  =  with 

So  if  we  take  a  =  ||(|w|  +  c)||oo  in  the  Lax-Friedrichs  splitting  (3.1), 


7=1  +  - 

1  TEfI^7e'(T) 

then  all  the  results  in  Section  3  will  hold.  The  following  CFL  condition  will  ensure  the 


positivity  of  w  +  2Ats(w). 


At  < 


P‘2 


if  uj  >  0 


2M2CJ  5 

- ,?}  ,  if  cu<0  • 

AMiuj  ’ 


(5.3) 


Consider  a  shock  tube  problem  for  the  reactive  flows  with  high  pressure  on  the  left  and 
low  pressure  on  the  right  initially  in  the  chemical  equilibrium  (uo  —  0).  The  initial  conditions 


are: 


(l Pl,  Tl)  =  (10001V/m2,  8000/1),  (pR,  Tn)  =  (1  N/m2,  8000/1 ) 


with  zero  velocity  everywhere  and  the  densities  satisfying 

Po  Po2  _  21  Pn2 

2  Ado  AIq-2  79  Ad]y 2 

where  Ado  —  0.016,  Mo2  =  0.032  and  Mat2  =  0.028.  The  initial  densities  of  O ,  O2  and  N2  are 
5.251896311257204  x  10“5,  3.748071704863518  x  10”5  and  2.962489471973072  x  10“4  on  the 
left,  and  8.341661837019181  x  10"8,  9.455418692098664 x  10”11  and  2.748909430004963 x  10“7 
on  the  right. 

High  order  schemes  without  the  positivity-preserving  limiter  may  blow  up  due  to  presence 
of  negative  pressure.  See  Figure  5.6  for  the  numerical  solution  of  the  WENO  scheme  at 
t  =  0.0001.  We  obtain  clean  and  grid  converged  solutions  for  this  test  case. 
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In  [10],  the  positive  real  number  to  avoid  denominators  becoming  zero  when  calculating 
nonlinear  weights  in  WENO  reconstruction  was  set  10~6.  For  this  particular  problem,  notice 
that  the  lowest  density  is  around  10-10,  so  we  adjust  the  positive  number  to  10~20.  Otherwise 
oscillations  may  emerge. 

6  Concluding  Remarks 

In  [20],  a  general  framework  was  established  to  construct  high  order  DG  and  finite  volume 
WENO  schemes  which  can  preserve  the  positivity  of  density  and  pressure  for  the  compressible 
Euler  equations  in  the  gas  dynamics.  In  this  paper,  we  have  shown  a  generalization  to  the 
finite  difference  WENO  scheme.  Extensions  to  multi-space  dimensions  are  straightforward. 

With  the  addition  of  the  positivity-preserving  limiter  in  this  paper,  which  involves  small 
additional  computational  cost,  to  the  finite  difference  scheme  (e.g.  ENO  and  WENO), 
the  numerical  solutions  will  always  have  positive  density  and  pressure  under  suitable  CFL 
condition.  We  have  tested  the  fifth  order  finite  difference  WENO  scheme  with  the  positivity¬ 
preserving  limiter  for  several  tough  problems.  Future  work  includes  the  generalization  to 
Navier-Stokes  equations. 
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(c)  Pn2 


(d)  Velocity 


(e)  Pressure 


Figure  5.6:  Example  5.6.  Fifth  order  WENO  scheme  with  positivity-preserving  limiter.  The 
solid  line  is  the  numerical  solution  of  Ax  =  The  symbols  denote  the  numerical  solution 

of  At  =  — — — 

or  lax  400Q. 
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